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ABSTRACT 


The  effort  reported  on  herein  has  concentrated  on  various  aspects  of  computing  unsteady 
flows  on  structured  dynamically  adaptive  meshes.  Various  wa>’s  of  coupling  the  mesh  motion 
to  the  flow  solution  were  studied.  Many  studies  were  performed  on  the  one-dimensional  un¬ 
steady  Euler  equations.  Without  the  complexity  of  multidimensional  phenomena  (both  phys¬ 
ical  and  numerical),  some  clear  trends  and  charaaeristics  were  identified  relating  to  accuracy 
of  the  computed  solutions  and  coupling  methodologies.  The  results  were  somewhat  surpris¬ 
ing.  Unsteady  flows  may  be  driven,  in  a  stationary  domain,  by  varying  a  boundary  condition  in 
time.  They  may  also  be  driven  by  varying  one  or  more  physical  boundary  positions  in  time, 
thus  simulating  a  moving  wall,  for  example.  This  work  has  looked  at  both  types  of  problems. 
The  moving  wall  problems  examined  were  both  one  and  two-dimensional  cases.  Issues  ex¬ 
amined  include  the  maximum  allowable  time  step  as  well  as  solution  accuracy  as  a  function  of 
amplitude  and  frequency  of  boundary  point  motion.  The  grid  point  positions  in  ail  cases  were 
obtained  as  a  solution  to  the  Euler-Lagrange  equations  resulting  from  a  variational  formula¬ 
tion.  The  grid  point  ^eeds  were  obtained  either  by  backward  differencing  the  grid  point  posi¬ 
tions  or  solving  the  grid  speed  equations  obtained  by  di^erentiating  the  grid  equations  tem¬ 
porally.  TWo  methods  of  solving  these  grid  speed  equations  with  differing  levels  of  difficulty 
were  employed  in  the  one-dimensional  problems  considered.  Tjvo-dimensional  problems 
studied  include  the  unsteady  flow  about  an  osdllating  wedge  body,  and  unsteady  flow  through 
a  CD  nozzle.  Some  preliminary  experience  in  computing  flows  with  a  shock  impinging  on  a 
thin  boundary  layer  region  is  also  presented. 
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I .  INTRODUCTION 


This  AFOSR/ARO  funded  effort  here  at  Iowa  State  University  started  on  November  1, 1989 
and  ended  on  October  31,  1992.  The  project  was  jointly  funded  by  AFOSR  and  ARO  and 
administered  by  AFOSR.  Early  discussions  with  both  AFOSR  and  ARO  scientists  resulted  in 
a  mutual  concurrence  to  diiect  the  research  effort  toward  the  eventuality  of  the  simulation  of 
a  moving  flexible  flight  vehicle/projectile.  The  unsteady  nature  of  both  the  fluid  flow  and  the 
geometry  in  this  case  necessitates  the  use  of.  and  a  thorough  understanding  of,  a  dynamically 
adapting  mesh.  The  mesh  for  such  a  problem  must  adapt  in  a  variety  of  ways.  It  must  adapt  to 
accommodate  boundary  surface  motion  as  well  as  to  provide  local  refinement  for  optimizing 
solution  accuracy.  Boundary  surface  motion  may  be  a  result  of  rigid  body  movement  of  all  or 
part  of  the  surface  and/or  it  may  be  a  result  of  local  deformation  of  surface  elements.  The 
significance  of  a  completely  time  accurate  algorithm  for  this  type  of  simulation  cannot  be 
overstated.  Issues  relating  to  the  feasibility  and  accuracy  of  this  type  of  simulation  are  the 
subject  of  this  effort  and  of  this  report. 

This  final  report  begins  with  a  review  of  the  project  goals.  Following  this  review,  a  summary  of 
effort  completed  is  presented  which  is  primarily  chronological.  Conclusions  are  presented 
which  represent  a  summary  of  all  significant  findings  as  they  are  expected  to  relate  to  the 
stated  direction  of  this  research.  Details  of  various  aspects  of  the  total  effort  are  presented  in 
the  appendix 
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II.  REVIEW  OF  CONTRACT  GOALS 

This  section  presents  a  list  summary  of  the  goals  defined  in  the  original  proposal  and  pres¬ 
ented  by  year  for  the  reader’s  convenience. 

A.  Yearl 

1)  Develop  animation  software  to  graphically  view  unsteady  motion  of  solution  and  gnd 

2)  Develop  dynamically  adaptive  mesh  scheme  for  the  wave  equation 

3)  Develop  mesh  scheme  for  the  invisdd  and  viscous  Burgers’  equation 

4)  Examine  errors  due  to  generalized  mappings 

5)  Examine  errors  associated  with  mesh  orthogonality 

6)  Examine  ability  to  track  shocks 

7)  Begin  development  of  unsteady  Euler  and  Navier-Stokes  codes 

B.  Year  2 

1)  Develop  mesh  scheme  for  one-dimensional  Euler  equations 

2)  Examine  errors  in  adapting  to  shock  tube  /  blast  wave  problems 

3)  Develop  mesh  scheme  for  two-dimensional  Euler  equations 

4)  Examine  ability  to  dynamically  adapt  to  two-dimensional  shocks 

C.  Year  3 

1)  Dcvelc^  mesh  scheme  for  two-dimensional  Navier-Stokes  equations 

2)  Implement  turbulence  model 

3)  Solve  flat  plate  flow  and  adapt  to  boundary  layer  gradients 

4)  Solve  shock-boundary  layer  impingement  problem 

5)  Solve  unsteady,  transonic  inlet  problem 
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Ill .  CHRONOLOGICAL  SUMMARY  OF  EFFORT 


A .  Graphics  Software 

The  first  task  was  to  develop  appropriate  simulation  software  for  visualization  of  unsteady 
results.  This  was  actually  an  ongoing  task  which  grew  in  complexity  as  the  need  to  visualize 
more  data  grew.  At  the  present  time,  the  simulation  software  runs  on  an  Apollo  workstation 
network  and  provides  graphical  output  of  grid  and  solution  evolution  for  the  1-D  and  2-D 
Euler  and  Navier-Stokes  equations.  The  graphical  output  is  in  the  form  of  a  computer  anima¬ 
tion  of  a  selected  sequence  of  frames  or  a  frame-at-a-time  view  of  a  user  chosen  variable. 

B  .  Adaptive  Mesh  On  1-D  Scalar  Problems 

A  grid  adaption  philosophy  was  presented  in  the  proposal.  Testing  of  this  philosophy  required 
test  cases  for  which  an  exact  solution  was  known.  This  required  that  the  dynamically  adaptive 
mesh  algorithm  be  implemented  on  some  simple  1-D  scalar  model  problems.  The  focus  here 
was  in  the  adaption  of  the  mesh  to  solution  changes  as  opposed  to  adaption  to  boundary  point 
motion.  The  adaption  procedure  chosen  was  one  based  on  a  variational  formulation  which 
was  an  extension  of  the  static  scheme  of  Brackbill  and  Saltzmanf  1  ]  into  the  dme  domain  by 
deriving  and  solving  a  companion  grid-speed  equation  simultaneously  with  the  unsteady 
model  equation.  A  similar  approach  based  on  a  set  of  Poisson  equations  was  presented  in 
12-4). 

The  need  to  begin  the  study  at  the  scalar  1-D  level  was  motivated  by  the  reality  that  in  order  to 
guarantee  success  in  the  multi-dimensional  problems  of  the  real  world,  it  is  essential  that  a 
complete  understanding  of  the  effect  of  introdudng  so  many  additional  degrees  of  freedom 
(doubled  for  the  scalar  1-D  problems  and  tripled  for  scalar  2-D  problems)  into  the  numerical 
problem  is  achieved.  The  effect  on  numerical  stability,  step  size  limitations,  accuracy,  wave 
formation,  wave  propagation,  etc.  were  of  interest.  Previous  work  has  mainly  adopted  a  “let’s 
see  if  we  can  do  it”  attitude  without  an  exhaustive  accuracy  assessment  based  on  known  exact 
solutions  and/or  experimental  data.  We  now  know  we  can  do  it.  What  we  do  not  know  is  how 
good  our  answer  is,  and  what  things  affect  this  answer,  and  in  what  ways.  In  fact,  in  general,  we 
do  not  even  know  how  to  define  “good”  in  the  present  context  for  problems  >»4iose  exact  solu¬ 
tions  are  not  known.  Consequently,  the  present  effort  sought  to  quantify  these  issues  to  the 
extent  possible. 

/ .  Numerical  Scheme 

A  number  of  test  cases  were  coded  and  run  for  the  wave  equation  and  Burgers’  equation 
which  included  moving  linear  and  non-linear  fronts  as  well  as  periodic  waves.  The  integra¬ 
tion  algorithm  used  was  the  standard  MacCormack  explicit  scheme.  This  was  chosen  as  a 
baseline  algorithm  due  to  its  popularity,  ease  of  extending  to  3-D,  and  ease  and  computation¬ 
al  efficiency  of  implementing  as  a  TVD  scheme  [  5  -  6  ]  (Note;  the  upwind  scheme  of  Warm- 
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ing  and  Beam  {  7  j  was  also  used  for  one  test).  1  he  ba.sic  algorithm  for  the  coupled  solver  is 
given  as  follows; 

1)  Given  u(x)  at  time  0,  and  an  initial  mesh  Xj 

2)  Apply  the  grid  law  at  time  0  to  generate  a  mesh  Xj  consistent  with  u^x) 

3)  Apply  the  gnd-speed  equation  at  time  0  to  generate  { x,  )j 

4)  Apply  the  predictor  step  of  the  algorithm  to  predict  u,  and  Xj  at  the  next  time,  t  +  At 

5)  Apply  the  grid  law  at  this  time  to  generate  a  mesh  Xj  consistent  with  this  Uj 

6)  Apply  the  grid-speed  equation  at  this  time  to  generate  ( x,  )j 

7)  Apply  the  a'>rrcctor  step  of  the  algorithm  to  improve  Uj  and  x,  at  this  time,  t  +  At 

8)  Apply  the  gnd  law  at  this  time  to  generate  a  mesh  x,  consistent  with  this  corrected  uj 

9)  Apply  the  grid-speed  equation  at  this  time  to  generate  ( x, ), 

10)  Go  to  step  4) 

A  modified  algorithm  was  also  tested  and  used  for  many  of  the  computations.  In  the  modified 
algorithm,  steps  5  and  8  are  eliminated  and  steps  6  and  9  are  replaced  by  the  line: 

6)  Apply  the  augmented  grid-speed  equation  at  this  time  to  generate  ( x,  )j 

where  the  word  “augmented”  is  used  to  indicate  the  use  of  the  grid  speed  control  law  ap¬ 
proach.  More  details  of  this  approach  are  presented  in  the  2-D  context  in  subsection  F  3  of 
this  section.  The  essence  of  it  is  that  the  grid  law  and  the  grid  speed  law  are  combined  into  one 
equation  with  feedback  control  characteristics  such  that  grids  obtained  by  integrating  grid 
speeds  maintain  there  form  as  solutions  of  the  grid  law  without  having  to  solve  the  grid  law 
explidtly.  This  saves  significant  CPU  time. 

Equations  of  the  form:  ^|udx+  j  h{x)*c  =  0  are  solved  in  a  finite-volume  context  on  a 

R  «R  R 

generalized  mesh  with  moving  control  volumes  using  a  IVD  version  of  MacCormack's 
scheme.  Here,  f=f(u)  and  h(x)  is  a  prescribed  forcing  function  which  is  sometimes  used  to 
make  the  problems  yield  solutions  which  tend  to  non  constant  functions  in  the  steady  state. 
The  grid  law  for  these  1-D  problems  has  the  form; 


(1) 


/'is  suggested  in  { 8  -  9  ],  the  dependent  variable  in  this  equation  may  be  the  arc  length  along 
the  function  w(x),  which  is  called  the  monitor  surface.  The  grid  speed  equation  corresponding 
to  this  grid  law  which  must  govern  the  proper  mesh  dynamics  is  given  by: 


iO 


<-2) 
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The  numerical  solutions  to  these  coupled  equations  for  the  test  problems  indicated  were  com¬ 
pared  directly  to  the  exact  solutions.  The  exact  error  provided  preci.se  information  as  to  the 
acceptability  of  the  various  weight  functions  and  corresponding  user  input  parameters. 

2 .  Coupling  Methods 

The  finite-volume  discrete  approximation  to  the  PDE  equation  has  grid  speed  in  it  a.s  part  of 
the  flux  evaluation.  The  grid  speed  equation  has  grid  speed  as  well  a.s  gradients  of  some  user 
selected  funaion  which  eventually  depends  on  gradients  of  the  PDE  solution.  These  two 
equations  are  inherently  coupled.  Three  basic  methods  of  .solving  these  equations  were  ex¬ 
amined  in  this  work.  They  are  called:  1)  differenced  speed  ( x, ).  2)  loosely  coupled  (differ¬ 
enced  Wt),  and  3)  strongly  coupled.  I  hese  methods  represent  drastically  different  degrees  of 
difficulty.  Method  1  is  by  far  the  simplest  and,  in  fact,  removes  the  need  of  even  having  a  grid 
speed  law.  The  grid  speeds  are  simply  determined  by  a  simple  backward  difference  of  the  grid 
point  positions.  Method  2  solves  the  grid  speed  law  but  the  quantity,  w^,  appearing  in  the 
equation  is  computed  from  a  backward  difference.  .Method  3  is  extremely  complex  in  that  the 
Wt  is  replaced  by  the  implied  sum  w^  =  w^  Uu  where  u,,  is  replaced  by  the  difference  approx¬ 
imation  to  the  governing  PDE  at  the  point  i.  This  renders  the  grid  .speed  equation  completely 
linear  (except  for  the  TVD  terms)  in  the  grid  speeds  and  allows  a  solution  for  the  grid  speeds 
with  the  PDE  approximation  and  all  of  the  associated  physics  of  the  subject  problem  to  influ¬ 
ence  the  results  in  a  time  accurate  fashion.  See  the  appendix  for  a  more  thorough  discussion 
of  the  coupling  procedures. 

3 .  Strong  Coupling 

The  strong  coupling  procedure  results  m  the  solution  of  a  nine-diagonal  system  of  equations. 
Substantial  effort  was  required  in  developing  this  system  and  finding  an  efficient  way  of  solv¬ 
ing  it.  The  system  character  is  a  result  of  using  a  smoothed  solution  to  drive  the  mesh.  This  is 
necessary  for  problems  which  amtain  abrupt  solution  features  such  as  shocks.  Two  smooth¬ 
ing  passes  result  in  a  nine-diagonal  system  whereas  three  passes  produce  an  eleven-diagonal 
system.  TVvo  passes  were  .sufficient  for  the  problems  studied. 

d .  Weight  Function 

The  mesh  adaption  scheme  u.sed  in  this  work  requires  a  minimum  of  two  user  input  parame¬ 
ters  which  control  the  amount  af  adaption  for  a  given  weight  function.  The  weight  function  is, 
however,  also  a  user  choice.  For  example,  the  user  may  elect  to  adapt  to  regions  where  the 
magnitude  of  the  gradient  of  a  chosen  function  is  large.  Or.  it  may  be  more  desirable  to  adapt 
to  regions  with  large  second  derivative,  or  perhaps  curvature,  or  a  combination  of  these 
things.  Such  a  generalized  weight  function  is  expressed  mathematically  a.s; 
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w  =  l+agi+bg2  +  cg3  +  dg4  4-  etc. 


(3) 


where  the  coefficients,  (a.b,c,cl,...)  are  user  prescribed  and  the  g,'s  are  squares  of  gradients, 
second  derivatives,  curvatures,  etc.  Even  after  this  decision  of  what  to  adapt  to  is  made,  the 
user  must  input  parameters  indicating  how  much  to  adapt.  While  some  work  has  been  done 
on  trying  to  automate  the  selection  of  the  input  parameters  for  certain  t)pes  of  weight  func¬ 
tions  [  10  ],  little  work  (with  notable  exceptions  (  8  -  9  j)  has  been  done  on  tr^ring  to  deter¬ 
mine  the  best  form  for  the  weight  function.  Because  of  this,  the  present  effort  was  initially 
concentrated  on  these  issues. 

5 .  Some  Accuracy  Tests 

The  error  in  the  numerical  solutions  at  a  point  is  defined  as  the  difference  between  the  com¬ 
puted  result  at  the  point  and  the  exact  result  at  the  point.  Another  measure  of  error  is  ob¬ 
tained  by  integrating  these  pointwise  errors  across  the  mesh  and  dividing  by  the  mesh  length. 
This  is,  of  course,  a  more  global  measure  of  error.  In  problems  with  shocks,  this  global  error 
can  be  totally  dominated  by  either  shock  speed  error  or  overshoot  error  when  TVD  schemes 
are  not  used.  Because  of  this,  care  in  the  interpretation  of  this  error  measure  must  be  exer¬ 
cised. 

Tests  for  these  1-D  model  equations  were  conducted  as  indicated  and  the  results  were  initially 
puzzling.  Results  indicated  the  following  general  conclusions. 

1)  TTie  adaptive  grid  solution  was  not  always  better  than  the  equivalent  fixed  gnd  solution. 

2)  The  adaptive  grid  solution  was  sometimes  unstable  when  the  foicd  grid  solution  was  not. 

3)  The  adaptive  grid  solution  could  always  be  made  better  thari  the  fixed  grid  solution  if  enough 
effort  was  spent  choosing  the  weight  function  and  the  adaption  intensity  parameters  but  the 
exact  solution  was  required  to  guide  these  user  choices. 

4)  The  user  choices  for  weight  function  and  adaption  intensity  parameters  which  provide  the  best 
result  for  one  problem  do  not.  in  general,  provide  the  best  result  for  a  different  problem. 

5)  The  user  choices  for  weight  function  and  adaption  intensity  parameter  which  provide  the  best 
result  for  one  numerical  algorithm  applied  to  a  given  problem  do  not  necessarily  provide  the 
best  result  for  a  different  algorithm  applied  to  the  same  problem. 

The  results  of  the  indicated  accuracy  tests  proved  that  simply  being  able  to  move  the  mesh 
points  dynamically  into  regions  with  prescribed  features  as  the  solution  changed  did  not  nec- 
e.ssarily  produce  a  numerical  solution  which  was  better  than  that  computed  on  a  fixed  mesh. 
This  meant  that  there  was  something  missing  in  the  understanding  of  what  was  actually  need¬ 
ed  to  constitute  a  “best  mesh”  or  a  “good  grid”  as  it  is  sometimes  referred  to.  Consequently, 
the  focus  of  the  effort  turned  to  this  issue. 

C  .  Errors  Due  To  Generalized  Mappings 
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The  following  path  of  reasoning  was  set  forward  to  address  this  issue:  Suppose  a  problem  of 
interest  requires  the  solution  of  a  partial  differential  equation  (PDE)  which  contains  a  deriva¬ 
tive  term,  f,,  where  f=  f(u(x))  with  u(x)  representing  the  dependent  variable  of  the  PDE.  The 
numerical  solution  of  such  a  problem  requires  first  a  discretization  of  the  x-domain  into  a 
sequence,  say,  {Xj,  i  =  0,1,...,!}  subject  to  the  constraint  that  xo<Xi  <  X2< ...  <X|,  and  second  a 
suitable  numerical  approximation  to  the  derivative,  f,.  The  objective  of  the  choice  of  the  par¬ 
ticular  sequence  {xj}  and  the  particular  discretization  used  in  approximating  f,  is  clearly  to 
minimize  the  error  between  the  numerical  value  of  f*  and  the  exact  value  of  f*  at  each  of  the 
mesh  points.  Clearly  then,  a  mesh  which  performs  this  function  is  a  “good  mesh”.  Moreover, 
a  mesh/discretization  combination  which  produces  zero  error  is  the  “best  mesh"  possible.  A 
numerical  procedure  was  then  developed  to  obtain  this  “best  mesh"  for  specified  functions 
and  derivative  approximations.  In  the  original  formulation,  the  actual  function,  f^x),  was  re¬ 
quired  for  the  numerical  procedure  suggested.  Details  of  this  are  found  in  the  appendix. 
Later  work,  however,  revealed  that  a  spline  fit  of  data  representing  the  function  was  adequate 
for  near  zero  error.  At  any  rate,  the  value  of  the  procedure  lies  not  specifically  in  its  practical 
applicability  (the  extent  of  which  remains  to  be  demonstrated),  but  rather  in  its  ability  to  pro¬ 
duce  categorically  the  “best  mesh”  for  any  differentiable  function  for  which  a  first  derivative 
expression  is  known.  This  does,  at  a  very  minimum,  provide  information  about  what  a  “good 
grid”  should  look  like  and  what  characteristics  it  should  have.  This  information,  however,  is 
valid  only  for  the  particular  discrete  approximation  used,  which  is.  in  most  cases  studied,  a 
central  difference.  It  is  possible  to  use  other  approximations  also  such  as  a  second  order  back¬ 
ward  difference  The  “best  mesh”  in  this  case  is,  strictly  speaking,  different  than  the  “best 
mesh”  in  the  central  difference  case.  The  important  thing  to  recognize  here  is  that  the  “best 
mesh”  depends  upon  both  the  nature  of  the  function  f(x)  and  the  choice  of  discretization  proce¬ 
dure  of  the  derivative,  f*.  Adaptive  mesh  procedures  tc  date  have  focused  on  only  the  first  of 
these.  In  theory,  one  could  determine  the  “best  mesh”  to  resolve  numerical  approximations 
to  other  types  of  terms  also,  such  as  the  quantity.  fx-pu„  found  in  the  viscous  Burgers'  equa¬ 
tion.  In  this  case,  however,  other  conditions  besides  no  inflection  points  must  likely  be  im¬ 
posed  to  guarantee  uniqueness  of  the  result,  if  such  uniqueness  is  even  possible. 

The  practical  usefulness  of  the  idea  suggested  here  is  unknown  for  multi-dimensional  multi- 
variable  problems.  Preliminary  work  shows  that  “best  grids”  may  be  obtained  in  2-D  for 
some  problems.  In  any  case,  the  lessons  learned  to  date  are  valuable  guides  in  helping  to 
formulate  new  weight  functions  and  new  adaption  schemes  which  incorporate  information 
regarding  both  the  nature  of  the  function  f(x)  and  the  discrete  approximation  used  to  represent 
its  derivatives. 

D  .  Extension  To  1-D  Unsteady  Euler  Problems 

The  basic  grid  and  grid-speed  laws  derived  for  the  scalar  model  problems  are  the  same  as 
those  for  the  1-D  Euler  equations.  Various  methods  of  coupling  the  mesh  dynamics  equa¬ 
tions  and  the  fluid  dynamics  equations  are  possible  as  mentioned.  These  methods  were  re- 
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ferred  to  as  methods  1, 2,  and  3.  Experiments  were  conducted  with  various  differencing  strat¬ 
egies  for  these  methods.  The  work  was  first  focused  on  1-D  shock  tube  computations  for  the 
problem  described  below.  In  these  computations,  the  weight  function  chosen  was  the  square 
of  the  density  gradient  plus  one.  This  caused  adaption  to  both  the  moving  shock  and  the  con¬ 
tact  discontinuity. 

1 .  Initial  Conditions 

The  initial  conditions  for  the  shock  tube  problem  were: 
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These  were  chosen  to  provide  for  a  comparison  with  the  results  in  (  6  ]  should  the  need  arise 
during  the  debugging  stages  of  code  development. 

2 .  Adaption  Errors  For  Shock  Tube 

In  order  to  continually  assess  the  quality  of  an  adaptive  mesh  solution,  certain  performance 
criteria  were  established.  The  absolute  pointwise  error  is  the  difference  between  the  com¬ 
puted  solution  at  a  mesh  point  and  the  exact  solution  at  that  point.  The  global  error  is  just  the 
integrated  pointwise  error  across  the  domain.  In  problems  with  shocks,  this  global  error  can 
be  totally  dominated  by  either  shock  speed  error  or  overshoot  error  when  TVD  schemes  are 
not  used  or  are  used  improperly.  Because  of  this,  care  in  the  interpretation  of  this  error  mea¬ 
sure  must  be  exercised, 

3 .  Results  for  Stationary  Boundaries 

The  findings  are  conclu.sive  in  some  respects  but  misleading  in  others,  for  the  cases  studied 
which  involved  only  solution  adaption  as  opposed  to  boundary  adaption.  For  problems  of  this 
type,  the  strongly  coupled  method  (method  3)  clearly  requires  the  most  computational  effort. 
The  differenced  speed  method  (method  1)  clearly  produces  the  most  accurate  results  of  the 
computations  performed  with  an  adapting  grid.  Figure  ( I )  shows  a  comparison  of  shock  tube 
solutions  for  this  case,  versus  the  exact  solution,  using  a  non-adapting  uniformly  spaced  grid, 
and  adapting  grids  using  methods  1,2,  and  3.  In  this  figure,  it  is  clear  that  method  1  produces  a 
solution  of  comparable  accuracy  (in  appearance)  to  the  solution  obtained  on  the  uniform  un¬ 
adapted  mesh.  It  also  appears  that  method  1  produces  better  results  for  this  ca«e  than  either 
of  methods  2  or  3.  It  was  of  more  than  casual  interest  that  the  solution  on  a  uniform  non¬ 
adapting  mesh  appeared  to  be  as  good  as  the  best  adaptive  grid  solution.  Indications  on  the 
shock  tube  problem  were  that  if  enough  points  were  used,  a  uniform  mesh  solution  was  more 
accurate  than  a  strongly  coupled  adaptive  mesh  solution.  It  appeared  that  there  was  some 
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critical  number  of  points  below  which  the  strongly  coupled  adaptive  mesh  solution  was  better 
than  the  uniform  mesh  solution.  For  the  shock  tube  problem  studied,  this  criticaJ  number  of 
points  was  around  80. 

Figures  (2)  and  (3)  support  this  notion.  Figure  (2)  illustrates  the  integrated  error  in  density 
versus  time  for  the  same  shock  tube  calculation  but  with  101  points  and  with  all  three  mesh 
speed  strategies.  The  one  labeled  “backward  P_tau"  refers  to  the  loo.sely  coupled  scheme  of 
method  2,  while  the  one  labeled  “backward  X  tau”  refers  to  the  differenced  speed  method 
(method  1).  It  is  apparent  from  this  figure  that  the  fully  coupled  method  produces  the  most 
error.  Interestingly,  this  method  also  takes  the  most  CPU  time.  The  smallest  error  and  small¬ 
est  CPU  time  occurs  for  the  differenced  speed  method.  Note  that  this  result  is  for  an  adaption 
constant  of  3.  In  these  figures,  the  error  resulting  from  no  adaption  is  presented  for  a  refer¬ 
ence,  The  minimum  error  in  both  of  these  figures  occurs  for  the  differenced  speed  method. 
As  pointed  out,  the  most  advantage  from  an  adaptive  mesh  scheme  should  come  from  compu¬ 
tations  on  a  coarse  mesh.  Figure  (3)  illustrates  the  same  computation  as  Fig.  (2)  except  51 
points  are  used  instead  of  101.  It  is  clear  from  these  figures  that  improved  solution  accuracy 
relative  to  the  uniform  mesh  accuracy  is  possible  with  adaption  on  a  coarse  mesh.  However, 
the  benefit  may  not  be  worth  the  price  in  CPU  time.  It  should  be  noted  that  a  typical  adaptive 
mesh  solution  took  0.0626  sec/step  as  compared  to  0.0 143  sec/step  for  the  uniform  mesh  case. 
In  addition,  due  to  the  step  size  restriction  imposed  by  the  finer  mesh,  it  took  2.4  times  as 
many  steps  to  reach  the  same  time.  The  results  here  make  it  clear  that  there  is  no  real  payoff 
for  this  problem  in  using  an  adaptive  mesh  based  on  a  weight  function  of  density  gradient. 
Better  weight  functions  undoubtedly  exist  and  substantially  more  improvement  is  possible  for 
other  types  of  problems.  However,  it  is  very  encouraging  that  the  adaptive  mesh  scheme  does 
an  excellent  job  of  accurately  tracking  shock  position  and  strength. 

It  would  seem  that  the  differenced  speed  method  is  the  method  of  choice  based  on  these  com¬ 
parisons.  As  mentioned  earlier,  these  results  are  valid  only  for  the  case  of  solution  adaption 
to  flows  with  shocks.  If  the  boundaries  of  the  domain  of  interest  were  set  into  motion,  these 
conclusions  need  not  apply.  This  possibility  was  tested. 

4.  Results  Jbr  Moving  Boundaries 

The  results  in  the  previous  section  were  bothersome.  Clearly,  mesh  adaption  with  ail  of  its 
overhead  did  not  improve  the  results  enough  to  justify  its  use.  The  previous  results  dealt  with 
solution  adaption.  The  boundaries  were  not  moving.  A  decision  was  made  to  study  the  effect 
of  setting  the  boundaries  in  motion.  Two  philosophies  for  doing  this  were  apparent.  An  oscil¬ 
lation  of  the  boundaries  could  be  interpreted  as  a  physical  disturbance  to  the  geometry  or  it 
could  be  viewed  as  simply  a  movement  of  the  computational  domain  boundaries  on  the  non¬ 
moving  geometry.  This  latter  scenario  was  chosen  for  the  test  case.  Each  boundary  of  the 
shock  tube  was  given  a  prescribed  sinusoidal  position  description.  For  the  left  end,  for  exam¬ 
ple,  the  left  boundary  point  moved  in  time  according  to: 
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xi{r)  =  M  sin(a)T  +  (p) 


(4) 


where  M,  io,  and  0  were  user  inputs.  The  right  boundary  was  treated  similarly.  Typically, 

M  was  set  to  some  constant  times  the  equivalent  uniform  Ax.  The  input  value  of  u)  was  chosen 
in  accordance  with  Tma*  by  requiring  a  selected  number  of  cycles  per  simulation  time.  The 
value  of  4)  was  0  for  all  cases  considered.  Many  computations  were  performed  with  method  1 
and  method  2  (several  hundred  cases).  From  the  data  generated,  the  integrated  error  in  den¬ 
sity  at  the  final  time  was  used  as  the  measure  of  accuracy  of  the  computed  solutions.  This 
integrated  error  is  plotted  versus  frequency  with  amplitude  as  a  parameter  in  Fig.  (4).  and 
versus  amplitude  with  frequency  as  a  parameter  in  Fig.  (5).  The  amplitude  is  M  in  these  fig¬ 
ures  and  the  frequency  is  the  number  of  cycles  in  the  fixed  amount  of  time  covered  by  the 
computations.  The  major  conclusion  from  these  figures  is  just  as  suspected.  The  benefit  of 
the  greater  coupling  between  the  grid  speed  equations  and  the  flow  equations  increases  as  the 
boundary  motion  increases.  The  error  present  at  the  lower  frequencies  and  magnitudes  is 
simply  a  result  of  the  truncation  error  of  the  numerics  used  to  approximate  the  grid  speed 
equations.  This  is  in  contrast  to  the  error  in  obtaining  the  grid  speeds  by  simple  backward 
difference.  As  the  integration  proceeds,  the  simple  backward  difference  of  method  1  allows  a 
build-up  of  error  which  could  be  interpreted  as  round-off.  The  higher  coupling  of  method  2 
inhibits  this  build-up  and  thus,  at  some  frequency  and  magnitude,  method  2  becomes  a  more 
accurate  approximation. 

This  completes  the  considerations  for  I-D  computations.  The  next  topic  is  the  extension  to 
2-D. 

E .  2-D  Considerations 

Other  considerations,  not  discussed  with  respect  to  the  1-D  examples  presented,  are  impor¬ 
tant  in  extending  to  multi-dimen.sional  problems.  One  of  the  most  obvious  is  the  issue  of 
orthogonality. 

/  .  Errors  Associated  With  Mesh  Non-OrthogonaiUy 

The  basic  concern  was  the  error  in  troduced  into  a  computation  as  a  result  of  non-orthogonal¬ 
ity  of  the  mesh.  In  a  dynamically  adapting  mesh,  this  concern  was  clearly  greater  since  the  real 
possibility  exists  that  cell  collapse  may  he  inadvertently  driven  by  the  adaption  scheme.  An 
attempt  was  made  to  study  the  way  in  which  non-orthogonality  influences  the  numerical  ap¬ 
proximation. 


a  .  Procedure  and  results 

The  basic  procedure,  outlined  here  and  detailed  in  the  appendix,  was  to  represent  the  deriva¬ 
tive  f*  in  a  2-D  generalized  coordinate  frame  with  derivatives  in  the  generalized  coordinate 
directions  approximated  by  central  differences.  The  f  quantities  were  then  expanded  about 
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the  reference  point  resulting  in  an  equation  of  the  form: 

=  Clfj,  +  C2fy  +  Cifxx  +  C^fyy  +  Csf^y  +  ... 

An  analysis  was  performed  to  see  how  the  grid  structure  locally  and  the  metric  evaluations 
influenced  the  Ci**.  In  summary,  the  basic  results  were  that: 

1.  metrics  should  be  evaluated  the  same  as  fj  ,  f,  to  zero  C2  and  make  C 1  unity. 

2.  minimizing  coordinate  line  curvature  and  stretching  also  mmimizes  numerators  of  C3,  C4,  C5. 

3.  minimizing  non-orthogonality  maximizes  denominators  of  C3.  C4,  C5. 

These  basic  results  are  expected  to  be  applicable  to  upwind  difference  approximations  as  well 
since  they  appear  to  be  common  sense  in  nature. 

F .  Dynamic  Adaption  For  2-D  Euler  Equations 

The  basic  grid  and  grid-speed  laws  for  2-D  problems  were  derived.  More  details  are  found  in 
the  appendix.  It  must  be  noted  that  the  equations  introduced  by  { 1  ]  are  dimensionally  incon¬ 
sistent  unless  dimensions  are  absorbed  into  the  user  input  parameters.  This  was  found  to  be 
inappropriate  resulting  in  difficulties  setting  the  values  of  the  parameters  for  a  desired  level  of 
grid  adaption.  The  parameter  values  for  a  given  level  of  adaption  were  problem  dependent. 
The  effectiveness  of  various  terms  in  the  performance  index  was  found  to  vary  over  the  do¬ 
main  in  problems  with  a  significant  clustering  requirement.  This  caused  the  effect  of  a  given 
term,  say  orthogonality,  to  be  more  significant  in  some  regions  than  it  was  in  others.  The  cure 
for  this  scaling  problem  was  to  redefine  the  performance  index  in  a  dimensionally  consistent 
manner  with  user  input  parameters  of  order  1.  The  Euler-Lagrange  equations  were  then 
re-derived  for  this  case  and  the  grid  speed  equations  were  derived  by  differentiating  these 
new  Euler-Lagrange  equations  with  respect  to  the  temporal  variable. 

7  .  Modified  Mesh  Law 

The  result  of  the  derivation  was  a  second  order  PDE  system  of  the  form. 

G(r)  =  Ar||  +  +  C?,^  +  Dr^  +  Er,,  =  0  ( 5 ) 

which  constitutes  the  grid  law.  The  boundary  conditions  are  mixed  Dirichlet  and  Neumann 
depending  upon  the  physical  conditions  required  of  the  mesh  at  the  boundaries. 

This  new  mesh  law  requires  specification  of  3  parameters,  X,,  Xq,  and  Xa  corresponding  to 
smoothness,  orthogonality,  and  adaption  terms.  Note  that  in  the  new  formulation,  these  may 
be  functions  of  time.  This  may  be  useful  for  problems  in  which  the  flow  structure  being 
adapted  to  diminishes  with  time  while  the  user  wishes  to  maintain  resolution  at  this  location 
during  the  decay  process.  These  parameters  are  nondimensional  and  have  the  value  of  unity  if 
orthogonality  and  solution  adaption  are  equal  in  importance  to  smoothness.  In  the  old 
scheme,  the  values  of  these  parameters  were  functions  of  the  problem,  the  number  of  mesh 
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points,  the  distribution  of  points,  etc.,  which  required  experimentation  with  various  values  in 
order  to  discover  an  acceptable  set  (if  one  even  existed). 

2 .  Corresponding  Mesh  Speed  Law 

The  result  of  the  temporal  differentiation  was  a  linear,  variable  coefficient,  second  order 
PDE  system  of  the  form: 


Gt(r)  =  +  C£,pj  +  D*z^  +  -f  f  *  =  0  ( 6  > 

which  constitutes  the  grid  speed  law.  In  this  equation,  z  =  Fj .  the  grid  speed  vector.  The 
boundary  conditions  are  mixed  Dirichlet  and  Neumann  for  this  equation  as  well,  depending 
upon  the  physical  conditions  required  of  the  mesh  speed  at  the  boundaries. 

3 .  Augmented  Mesh  Speed  Law 

The  mesh  law  and  mesh  speed  law  were  combined  into  one  equation  called  the  “augmented” 
mesh  speed  law.  This  equation  is  written  symbolically  as: 

G^r)  +  .IcGCr)  =  0  (7) 

where  Xc  is  a  user  spedhed  feedback  constant  which  was  typically  set  to  a  number  between  50 
and  500.  This  number  is  of  order  1/At.  This  equation  is  used  as  a  feedback  control  law  which 
removes  the  need  for  solving  the  grid  law.  When  the  grid  speeds  obtained  from  solving  this 
“augmented”  grid  speed  equation  are  simply  integrated  in  time  to  obtain  the  grid  at  the  new 
time,  this  new  grid  automatically  satisfies  the  grid  law. 

4 .  Solution  Strategies 

Based  on  the  findings  in  1-D  regarding  the  strongly  coupled  procedure  and  based  on  the 
enormous  complexity  and  CPU  requirement  of  this  method,  it  was  rejected  as  a  candidate  for 
the  2-D  computations.  The  two  methods  called  methods  1  and  2  were  employed.  The  TVD 
version  of  MacCormack’s  method  was  extended  to  2-D  and  used  as  the  numerical  integrator. 

5 .  Oscillating  Wedge 

The  first  2-D  test  on  Euler  equations  is  an  oscillating  wedge/plate  interaction  problem  as 
shown. 

The  8(t)  is  a  user  prescribed  function  as  is  the  point  about  which  rotation  is  to  take  place.  The 
resulting  shock  pattern  will  normally  oscillate  between  a  regular  reflection  pattern  and  a 
Mach  reflection  pattern.  The  incident  shock  may  be  planar,  curved,  or  non-existent  (it  could 
be  an  expansion)  depending  upon  the  magnitude  of  the  oscillation.  The  triple  point  solution 
for  the  Mach  reflection  problem  is  known  for  steady  problems.  The  regular  reflection  solu¬ 
tion  for  the  planar  shock  is  also  known  for  the  steady  problem.  This  information  helped  to 
validate  some  of  the  unsteady  solutions. 
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The  test  run  involved  a  general  8(t)  as  shown  below. 


The  high  and  low  values  of  8  were  input  as  well  as  the  duration  times  of  each  phase  of  the 
curve.  This  defines  a  general  periodic  input  which  can  be  a  sinusoid  as  well  as  a  more  general 
function.  An  example  of  a  simple  compression-expansion  process  follows. 

Results  were  computed  but  are  not  shown  here  for  a  20  degree  wedge  in  a  Mach  2  flow  with 
two  schedules  back  to  back  which  cause  the  wedge  to  start  at  8|  =  -10  degrees  which  corre¬ 
sponds  to  no  flow  deflection  on  the  lower  side.  The  wedge  then  moves  along  a  sine  curve  to 
8h  -  Othusproducinga  maximum  compression  deflection  of  10  degrees.  This  is  followed  by  a 
retreat  back  to  8  =  -10  degrees  and  then  on  to  an  expansion.  The  results  of  a  computation  of 
this  sort  involve  over  a  thousand  frames  of  data,  many  of  which  are  viewed  on  the  workstation 
as  (in  some  cases,  after)  the  solutions  are  computed.  The  management  of  this  much  data  is  an 
enormous  task.  Graphical  procedures  have  been  developed  to  enable  the  management  and 
processing  of  this  data.  The  data  can  now  be  made  into  an  animation  video  for  use  in  presen¬ 
tations. 
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One  of  the  cases  which  was  run  and  is  presented  here  started  the  lower  wedge  surface  for  the 
20  degree  wedge  parallel  to  the  freestream  and  then  the  body  was  plunged  downward  until  a 
prescribed  compression  angle  on  the  lower  surface  was  achieved.  This  case  is  similar  to  the 
previous  one  except  that  when  the  compression  angle  was  achieved,  the  body  remained  at  that 
angle.  The  unsteady-  flowfield  in  this  case  started  as  a  parallel  supersonic  freestream  and  be¬ 
gan  a  compression  process  which,  after  much  shock  development  and  interaction,  ultimately 
caused  the  flow  to  reverse,  first  at  the  downstream  boundary,  and  then  proceed  to  the  left. 
This  process  continued  until  the  entire  flow  was  reversed.  This  is,  in  effect,  a  simulation  of  the 
choking  process  in  a  supersonic  flow.  The  animation  of  this  flow  evolution  was  indeed  fasci¬ 
nating.  A  VHS  format  video  of  this  animation  was  made  and  a  copy  of  it  is  available  upon 
request.  It  is  not  in  color  as  the  hardware  available  at  the  time  prevented  this  possibility.  Two 
pages  of  selected  frames  illustrating  the  flow  evolution  are  presented  in  Figs.  (6)  and  (7).  The 
first  of  these  illustrates  Mach  contours  while  the  second  shows  velocity  vectors.  Progress  is 
from  left  to  right. 


In  the  computation  of  this  case,  it  became  apparent  that  the  explicit  CFt .  condition  alone  was 
not  always  sufficient  for  the  determination  of  the  stable  time  step.  Some  effort  was  spent 
studying  this  in  more  detail.  A  list  of  the  relevant  time  scales  (on  the  mesh  scale)  was  com¬ 
piled-  The  two  conventional  ones  which  come  from  an  Eigenanalysis,  are  given  by: 


Ar, 


Ax 

|u-Xr|  +  a 


Afy 


Ay 

v-yrl  +  a 


(8) 


In  these  equations,  (u,v)  and  (xT,yT.)  represent  the  components  of  the  velocity  and  the  ^ed  of 
the  mesh  point  in  the  x,y  directions  respectively,  and  ’a’  represents  the  acoustic  velocity.  A 
third  time  scale  which  seems  intuitive  for  this  problem  is  given  by 


Ay  _  Ay 
(Vbodyl  ■  R\(0\ 


where  («)  is  the  angular  velocity  of  the  body  motion  and  R  is  the  distance  from  the  pivot  point  to 
the  point  in  question.  This  time  scale  is  distinctly  different  than  Aty  above  since  v-yT  vanishes 
at  the  body.  This  time  scale  results  from  the  restriction  that  the  body  motion  can  not  cause  a 
mesh  point  to  sweep  aaoss  more  than  one  existing  mesh  cell  in  a  single  time  step.  This  restric¬ 
tion  is  consistent  with  the  notion  of  time  accurate  unsteady  flow  simulation.  It  appears  that  it 
is  the  relationship  of  this  time  scale  to  the  smallest  of  the  other  two  which  determines  if  there 
is  any  advantage  to  using  implicit  methods.  If  At^  is  significantly  larger  than  the  minimum  of 
the  other  time  scales,  then  implicit  methods  are  likely  to  be  of  significant  benefit.  A  more 
general  statement  of  this  idea  is: 
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and  the  minimum  of  these  stepsizes  determines  the  significant  timescale  to  be  compared  with 
the  conventional  ones  from  eigenanalysis. 

6.  CD  nozzle 

A2-D  computation  of  flow  through  a  converging-diverging  nozzle  was  made  and  compared 
with  a  solution  of  the  model  problem  of  quasi  1-D  flow  through  a  variable  area  duct.  The  case 
of  interest  is  one  with  a  combination  of  inlet/exit  conditions  which  cause  a  shock  to  form  just 
downstream  of  the  throat.  The  exit  flow  Mach  number  is  subsonic  in  this  case  which  requires 
specification  of  one  flow  condition  there.  The  variable  chosen  was  static  pressure.  This  exit 
pressure  was  varied  sinusoidally  thus  causing  the  shock  position  to  oscillate  in  the  nozzle.  An 
interesting  aspect  of  this  problem  was  that  if  the  back-pressure  oscillation  was  strong  enough 
and  rapid  enough,  it  appeared  that  the  pressure  wave  generated  by  this  boundary  condition 
was  nearly  a  shock  which  traveled  upstream  and  interacted  with  the  shock  already  present  in 
the  nozzle. 

G  .  Dynamic  Adaption  For  2-D  Navier-Stokes  Equations 

Extension  of  the  dynamic  adaption  procedures  to  viscous  flows  was  performed.  Some  simple 
test  cases  were  computed.  These  included  i  flat  plate  boundary  layer  calculation  in  both  sub¬ 
sonic  and  supersonic  flow,  and  a  shock-boundary  layer  interaction  problem  reported  on  by 
Hakkinen  et  al.  (  11  ].  Results  for  this  latter  case  are  presented  here. 

/  .  Considerations  for  Viscous  Flows 

Dynamic  adaption  for  viscous  flow^  was  and  remains  a  real  challenge.  Grid  points  are  re¬ 
quired  to  cluster  near  the  wall  in  boundary  layer  regions  while  at  the  same  time  clustering  to 
features  in  the  outer  flow  such  as  shocks.  It  was  unclear  what  weight  function  would  provide 
the  best  of  both  worlds.  One  ba.sed  on  Mach  number  gradient  was  used. 

Further  complications  arose  due  to  the  inadequacy  of  simple  turbulence  models  to  correctly 
model  flow  physics  in  cases  with  turbulence.  This  ever-present  reality  was  exacerbated  by  the 
additional  degrees  of  freedom  present  in  the  system  of  equations  due  to  the  dynamically 
adapting  grid,  and  by  the  difficulty  in  getting  the  grid  to  adapt  enough  in  the  boundary  layer 
region  without  over-adapting  at  the  shock,  even  in  laminar  boundary  layer  problems  where 
wall  clustering  requirements  are  significantly  less  than  those  for  turbulent  boundary  layers.. 

2 .  Shock-Boundary  Layer  Interaction 

The  supersonic  flow  over  a  flat  plate  with  an  impinging  shock  was  studied  in  Ref.  [11].  The 
flow  conditions  for  this  test  case  were:  M<»  =  2,  Reo©  =  2.96(10)^  Too  =  250‘’K,  and  L=  1  me- 
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ter.  The  ■'omputations  for  this  case  were  made  on  both  a  clustered  static  mesh  and  a  dynami¬ 
cally  adaptive  mesh.  The  static  mesh  and  the  converged  dynamically  adaptive  mesh  are  shown 
in  Fig.  (8).  It  is  clear  that  there  is  significantly  less  clustering  in  the  dynamically  adapting  mesh 
in  the  boundary  layer  region.  Figure  (9)  illustrates  the  pressure  contours  for  both  the  static 
mesh  and  the  dynamically  adaptive  mesh.  The  shock  resolution  is  clearly  improved  by  the 
adaptive  mesh.  However,  surface  pressure,  for  example,  shown  in  Fig.  (10),  demonstrates 
that  the  adaptive  mesh  solution  misses  the  separation  point  and  the  experimental  data,  in  gen¬ 
eral,  more  than  the  static  mesh  solution.  It  is  felt  that  this  discrepancy  is  mainly  due  to  lack  of 
resolution  of  the  boundary  layer  region  by  the  adapting  mesh. 
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IV.  CONCLUSIONS 


Dynamic  adaption  of  a  mesh  is  a  result  of  either  boundary  motion  or  solution  adaption  or 
both.  Solution  adaption  may  be  a  result  or  boundary  motion,  boundary  condition  unsteadi¬ 
ness,  or  just  inherent  solution  unsteadiness.  The  present  research  results  suggest  that  un¬ 
steadiness  resulting  from  boundary  motion  may  require  a  greater  level  of  coupling  between 
grid  motion  prediction  and  flow  solution  prediction  than  unsteadiness  resulting  from  simple 
moving  waves.  It  also  is  very  clear  that  dynamically  adaptive  meshes  are  very  CPU  intensive 
and  should  be  used  only  for  problems  on  which  there  use  is  needed.  Information  is  presented 
in  this  report  which  guides  the  researcher  regarding  when  adaptive  meshes  are  useful  and 
when  they  just  produce  exr  "'ss  baggage.  Information  is  also  presented  which  guides  the  choice 
of  weight  function  in  the  context  of  the  definition  presented  for  constructing  the  “best  mesh,” 

This  work  tested  some  ideas  regarding  the  significance  of  various  time  scales  present  in  a  flow 
problem  with  dynamic  motion.  The  results  from  this  testing  suggest  a  means  of  guessing  when 
it  is  benefidai  to  use  implicit  methods  instead  of  explidt  ones. 

The  various  methods  of  coupling  the  mesh  dynamics  and  the  flow  dynamics  is,  of  course,  a 
crudal  issue.  The  strong  coupling  procedure  is  the  most  elegant,  and  the  most  rigorous.  It 
does  not  produce  the  best  results  in  the  cases  discussed  in  this  report,  but  as  mentioned,  is 
ejq)ected  to  be  superior  for  problems  in  which  there  is  complex  and/or  rapid  mesh  motion  due 
to  complex  flow  feature  evolution  or  complex  boundary  motion.  Unfortunately,  it  is  unques¬ 
tionably  the  most  expensive.  Perhaps  the  greatest  disadvantage  of  the  strong  coupling  proce¬ 
dure  is  that  changing  the  nature  of  the  weight  function  requires  a  significant  analysis  effort 
and  a  non-trivial  code  modification  effort.  The  loosely  coupled  procedure  (method  2)  does 
not  suffer  from  this  difficulty  and  consequently  is  the  method  of  choice  for  problems  requiring 
dynamic  adaption. 

The  mesh  adaption  procedures  explored  in  this  research  are  ideally  suited  for  use  in  problems 
with  moving  as  well  as  deforming  boundaries. 

The  work  has  been  carefully  guided  to  produce  numerical  procedures  which  are  ail  extend¬ 
able  to  3-D.  This  extension  is  straightforward. 
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Figure  (1).  uniform,  method  1.  2.  and  3,  and  exact  solutions. 
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Figure  (2).  integrated  density  error  versus  time  (101  points) 
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Figure  (3).  integrated  density  error  versus  time  (51  points) 


»  BtfTMk^O  O 

»  • 

t  2 

.  1  • 
I  O 

>  pU.m^gi  O  € 
.  2 

r  8 


18  47-Ot 


frequenqf,  cycles 


Figure  (4).  integrated  density  error  versus  oscillation  frequency  for  methods  1  and  2 
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Figure  (5).  integrated  density  error  versus  oscillation  magnitude  for  methods  1  and  2 
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Figure  (6).  Mach  contour  history  for  oscillating  wedge  flow 
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Figure  (7).  velocity  vector  history  for  oscillating  wedge  flow 
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Figure  (8).  static(a)  and  dynamically  adapted(b)  meshes  for  shock-boundary 

layer  problem 
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Figure  (9).  pressure  contours  for  static(a)  and  dynamically  adapted(b) 
meshes  for  shock-boundary  layer  problem 
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Figure  (10).  plate  pressure  cx)efficient  for  shock-boundary  layer  interaction 
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VIII.  APPENDIX 


A  .  Adaptive  Mesh  On  1-D  Scalar  Problems 

The  finite-volume  form  of  the  governing  equation  as  applied  to  a  moving  control  volume  has 
terms  involving  the  grid  speed.  The  solution  at  the  next  time  step  requires  the  values  of  the 
grid  speed  at  each  of  the  mesh  points.  The  grid  speed  equation,  Eq.  ( 2 ),  has  the  quantity,  w^, 
in  it,  in  addition  to  the  primary  unknown,  Xt-  The  w  in  this  equation  depends  on,  say  gj  (see 
Eq.  (  3  )),  which  might  be  something  like  (u^^  for  example.  Thus,  the  t  derivative  of  w  ulti¬ 
mately  depends  on  u^,  which  in  turn  depends  on,  among  other  things,  Xr.  The  point  is  that  this 
coupling  between  the  governing  equation  and  the  grid  speed  equation  may  be  done  numeri¬ 
cally  in  more  than  one  way.  The  next  section  details  the  two  methods  used  in  this  work. 

1 .  Methods  of  CoupUng 

a  .  Looseh  Coupled 

In  this  method,  Eq.  ( 2  )  was  discretized  by  using  central  differences  for  all  .spatial  discretiza¬ 
tion  and  first  order  backward  differences  for  w^.  The  discrete  system  of  equations  was  then 
solved  by  a  point  relaxation  scheme.  Line  relaxation  was  used  on  some  tests  but  did  not  con¬ 
verge  as  well  as  the  point  schemes  for  several  problems.  The  point  schemes  were  of  Jacobi 
and  Gauss-Siedel  types,  as  well  as  a  Newton  iteration  scheme.  The  point  Newton  iteration 
scheme  performed  the  best  overall. 

b.  Strondv  Coupled 

This  method  is  quite  complex  and  required  extensive  coding,  even  for  the  1-D  case.  It  is  the 
most  eloquent  method  of  coupling,  however.  The  essence  of  it  is  that  the  dependencies  indi¬ 
cated  in  the  open  paragraph  of  VIII  A  are  followed  in  a  most  general  way.  This  is  done  by 
defining  a  function,  P,  as  P  =  wg/w.  This  makes  the  grid  law  appear  simpler  as: 

+  x^P  =  0 .  ( 10 ) 

The  corresponding  grid  speed  law  is; 

+  (Xr)^P  -b  X^P,  =  0.  (ll) 

Pj  represents  P  at  the  point  i  and  is  expressed  as  a  function  of  x  and  u  .  These  vectors  repre¬ 
sent  all  ofthe  x  and  u  values  in  the  mesh.  Thus,  (P,),  =  (P,)iiu,  +  (Pi>^.  The  notation  on  the 
right  hand  side  of  this  equation  implies  a  summation  ofthe  terms  corresponding  to  each  ofthe 
elements  of  the  x  and  u  vectors.  Of  course,  the  Ur  is  replaced  with  terms  involving  x,  from 
the  governing  equation  for  the  problem.  The  matrices  (Pi)j  and  (P,>f  are  a  function  of  the 

exact  nature  of  the  chosen  weight  function,  w.  To  complicate  matters  ftjrther,  the  actual  u 
used  in  this  equation  was  a  smoothed  u  rather  than  the  computed  one.  This  was  necessary 
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since  the  computed  u  often  had  anomalies  due  to  shocks,  etc.,  which  needed  to  be  washed  out 
before  using  it  to  drive  the  mesh.  This  method  produces  a  grid  speed  equation  which  is  linear 
in  the  unknown  <  except  for  the  TVD  terms.  This  equation  was  discretized  with  central 
differences  for  all  spatial  derivatives  except  those  resulting  from  the  substitution  for  Ut. 
Those  were  treated  identical  to  their  treatment  in  the  di.scretization  of  the  governing  equa¬ 
tion.  The  net  result  of  this  procedure  was  a  linear  sv-stem  of  equations  to  solve.  This  was  a 
nine-diagonal  system  of  equations  for  which  a  special  solver  was  written.  The  bandwidth  of 
this  system  was  a  function  of  how  many  smoothing  passes  were  used.  No  smoothing  resulted 
in  a  pentadiagonal  system.  The  nine-diagonal  was  a  result  or  two  smoothing  passes.  This 
system  of  equations  was  also  solved  by  various  point  relaxation  .schemes  with  varying  degrees 
of  success. 

B  .  Errors  Due  To  Generalized  Mappings  (1-D) 

This  section  expands  section  III  C  and  defines  the  meaning  of  a  “good  grid".  Suppose  a  prob¬ 
lem  of  interest  requires  the  solution  of  a  partial  differential  equation  (PDE)  which  contains  a 
derivative  term,  f*,  where  f  =  f(u(x))  with  u(x)  representing  the  dependent  variable  of  the  PDE. 
The  numerical  solution  of  such  a  problem  requires  first  a  discretization  of  the  x-domain  into  a 

sequence,  say,  (xj,  i  =  0,1 1}  subject  to  the  constraint  that  xq  <  xi  <  X2  < ...  <  xi,  and  second  a 

suitable  numerical  approximation  to  the  derivative,  f,.  The  objective  of  the  choice  of  the  par¬ 
ticular  sequence  {x,}  and  the  particular  discretization  u.sed  in  approximating  f,  is  clearly  to 
minimize  the  error  between  the  numerical  value  of  f*  and  the  exact  value  of  f,  at  each  of  the 
mesh  points.  Clearly  then,  a  mesh  which  performs  this  function  is  a  “good  mesh".  Moreover, 
a  mesh/discretization  combination  which  produces  zero  error  is  the  “best  mesh”  possible. 

/  .  The  Best  Mesh  (A  concept) 

Since  the  same  PDE  may  yield  many  different  solutions  depending  on  the  particular  initial 
and  boundary  conditions  imposed,  the  optimal  sequence  {x,}  (i.e.,  the  best  mesh)  generally 
will  be  different  for  different  problem‘;  In  addition,  the  discrete  approximation  to  f,  could 
also  be  different  for  different  problem.^.  This  is  a  common  occurrence  in  many  currently  pop¬ 
ular  schemes  in  which  the  discretization  is  determined  based  on  local  information  in  the  solu¬ 
tion.  Upwind  schemes  and  schemes  ba.sed  on  a  solution  s  total  variation  (TV)  (e.g..  TVD 
schemes)  are  examples  of  this.  Prevailing  thought,  therefore,  is  usually  to  first  decide  upon  a 
discretization  procedure  for  the  PDE  of  interest.  Then  a  procedure  is  developed  for  adjusting 
the  mesh  point  positions  to  improve  the  resolution  of  the  computed  solution.  These  are  gen¬ 
erally  treated  as  separate,  if  not  unrelated,  problems.  Mesh  point  adjustment  schemes  are 
usually  based,  in  part,  on  a  user's  intuition.  For  example,  it  is  common  to  find  attempts  at 
concentrating  mesh  points  in  regions  with  large  values  of  jf,!  independent  of  the  particular 
discretization  scheme  used.  Such  intuitively  driven  mesh  adjustment  schemes  work  some¬ 
times  for  some  f{x)  distributions  but  not  for  others.  The  following  analysis  helps  to  explain 
why. 


35 


Given  a  differentiable  function  f{x),  consider  the  following  finite-difference  approximation 
to  the  first  derivative  of  this  function  at  the  point  x,; 


fx(X;) 


f|  ->-  1  ~  f|-l 


f  12  ) 


where  fi+ 1  means  f(x,  +  f).  etc.  This  approximate  expression  may  be  written  as  a  stria  equality 
provided  the  Xj  on  the  left  hand  side  is  interpreted  as  the  value  of  x  between  Xj.]  and  X|4- 1  at 
which  Eq.  (  12  )  becomes  exact  in  the  sense  of  the  mean  value  theorem.  Thus 


U\) 


fl-*- 1  ~  fj-l 


(  13) 


The  problem  shaping  up  here  is  to  find  a  point.  Xj,  between  each  pair  of  bounding  points  Xj.i 
and  Xj  + 1  which  causes  the  finite-difference  approximation  at  the  point  Xj  to  produce  the  exaa 
funaion  derivative  at  this  point.  Note  that,  strialy  speaking,  the  exact  function  derivative 
expression  is  required  for  such  an  approach  to  proceed.  Further  note  that  even  though  the 
determination  of  Xj  requires  an  iterative  solution  to  a  non-linear  equation  in  general,  a 
unique  result  is  guaranteed  by  the  mean  value  theorem  provided  the  sequence  {fj}  on  the 
mesh  {xj}  adequately  represents  the  true  function  nature.  More  specifically,  fjx  must  be  of 
one  sign  between  each  pair  of  bounding  points  x,.j  and  Xi+ 1.  This  requirement  presents  a 
problem  for  intervals  containing  inflection  points  e,xcept  that  any  given  curve  can  be  broken 
into  a  sequence  of  curves,  none  of  v.  hich  have  inflection  points  on  their  interior.  The  break¬ 
points  of  this  .sequence  are,  of  course,  the  inflection  points  themselves. 

Since  the  actual  function,  f,(x),  is  required  for  the  numerical  procedure  suggested  above,  it 
would  seem  a  fiaiitless  procedure  tor  prohlemsof  practical  interest.  However,  the  value  ofthe 
procedure  lies  not  specifically  in  its  practical  applicability  (the  extent  of  which  remains  to  be 
demonstrated),  but  rather  in  its  ability  to  produce  categorically  the  “best  mesh”  for  any  differ¬ 
entiable  function  for  which  a  first  derivative  expression  is  known.  This  will,  at  a  very  mini¬ 
mum,  provide  information  about  what  a  “good  grid"  should  look  like  and  what  charaaeristics 
it  should  have.  This  information,  however,  will  be  valid  only  for  the  particular  discrete  ap¬ 
proximation  used  in  Eq.  ( 13  ),  which  is,  in  this  case,  a  central  difference.  It  is  possible  to  use 
other  approximations  also.  For  example,  a  second  order  backward  difference  would  require  a 
solution  of  the  equation; 
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fx(x,) 


-H  f,.2 
3X|-4Xj_;  +  Xi-2 


(  14  ) 


for  the  “best”  location  of  the  point  x,.  The  “best  mesh"  in  this  case  is  different  than  the  "best 
mesh"  in  the  central  difference  case.  The  important  thing  to  recognize  here  is  that  the  "best 
mesh”  depends  upon  both  the  naiure  of  the  function  f(x)  and  the  choice  of  discretization  proce¬ 
dure  of  the  derivative,  f*.  Adaptive  mesh  procedures  to  date  have  focused  on  only  the  first  of 
these. 

2  ,  Some  First  Derivative  Examples 

The  actual  iterative  procedure  used  for  solving  Eqs.  (  13  )  and  (  14  )  is  now  outlined.  Differ¬ 
ent  types  of  functions  are  used  to  illustrate  the  ability  of  the  method  to  produce  the  “best 
mesh”  for  both  central  difference  approximations  and  one-sided  difference  approximations. 
As  a  beginning,  consider  the  iterative  form  of  Eq.  (  13  ): 


This  form  implies  a  double  iteration  nested  with  the  local  interval  iteration  inside  of  an  outer 
iteration  spanning  the  mesh  from  lowest  index  to  highest  index.  This  iterative  procedure  is 
equivalent  to  a  Newton  iteration  embedded  within  a  cyclic  mesh  sweep  outer  iteration.  For 
each  value  of  the  mesh  index,  i,  and  starting  from  an  initial  condition  for  the  sequence  on  i, 
{Xj}*' (k  =  0),  generate  as  the  limit  of  the  convergent  sequence  on  m,  {xf’’^ ^ .  This 

iterative  process,  of  course,  requires  f(x)  to  be  of  class  but  only  because  of  the  choice  to 
embed  a  Newton  iteration  in  the  procedure.  An  equivalent  method  can  be  applied  to 
Eq.  (  14  )  for  a  one-sided  difference  approximation.  It  should  be  observed  that  the  proce¬ 
dure  must  fail  if  the  function  f(x)  is  a  straight  line,  since  in  that  case,  it  no  longer  matters  where 
the  mesh  points  are  placed  from  the  point  of  view  of  evaluating  a  first  derivative  with  a  finite- 
difference.  The  resulting  difference  approximation  will  always  produce  the  exact  derivative 
regardless  of  where  the  mesh  points  are. 

Selected  examples  of  the  application  of  this  procedure  for  various  functions  are  now  pres¬ 
ented.  The  examples  are  selected  to  demonstrate  some  significant  results. 

a  ■  ffx)  =  sinfTLxi 

The  function  sinf-nx)  serves  as  the  first  test  case.  Figure  ( 1 )  shows  the  function  along  with  the 
two  “best  mesh”  distributions.  Note  that  the  points  tend  to  cluster  toward  the  peak  region  for 
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both  finite-difference  approximations.  Further  note  that  the  two  finite-difference  approxi¬ 
mations  produce  different  “best”  meshes  and  that  the  mesh  for  the  one-sided  scheme  is  direc¬ 
tionally  biased  as  one  would  expect.  A  traditional  adaptive  mesh  approach  would  adapt  to 
some  weight  function  chosen  by  the  user.  This  might  be  function  gradient,  second  derivative, 
curvature,  or  some  combination  of  these.  These  three  functions  are  plotted  with  f{x)  in  Fig. 
(2).  Note  that  a  weight  function  based  on  either  f.  f„.  or  curvature  would  at  least  provide  the 
proper  clustering  direction  for  this  case  whereas  one  based  on  f,  would  not.  The  determina¬ 
tion  of  the  “best  mesh”,  however,  requires  more  than  just  directional  information. 

b ,  /Crj  =  fcrif  QV  t  dii 

In  this  example,  a  =  .8(p  4- 4),  b  =  -1.6(p  +  4).  c  =  p.  d  =  -.2(p-16),  where  p  is  a  user  input  pa¬ 
rameter  which  is  equal  to  half  of  the  value  of  f,,  at  x  =  0  and  controls  the  nature  of  the  quartic. 
The  first  case  is  for  a  user  input  of  p  =  -8.  The  function  along  with  the  two  “best  mesh”  distri¬ 
butions  are  shown  in  Fig.  (3).  This  figure  clearly  illustrates  that  the  points  in  the  “best  mesh" 
for  this  problem  tend  to  cluster  toward  the  boundaries.  Although  this  function  exhibits  some 
of  the  same  qualitative  characteristics  as  the  sin(iTx)  function  as  seen  in  Fig.  (4),  the  points  in 
the  “best  mesh”  do  the  exact  opposite  of  what  they  do  in  the  quartic  case.  Obviously,  an  adap¬ 
tive  mesh  procedure  which  clustered  toward,  say,  high  first  derivative  regions  would  perhaps 
improve  the  resolution  of  the  derivative  f*  for  the  quartic  function  but  would  cause  poorer 
resolution  for  sin(Trx). 

The  final  demonstration  is  made  with  p  =  -4.  Figure  (5)  shows  the  “best  mesh”  for  this  case 
while  the  possible  weight  functions  are  shown  with  f(x)  in  Fig.  (6).  Note  that  for  this  case, 
there  is  no  clustering  at  all  for  either  of  the  finite-difference  approximations.  It  turns  out  that 
this  case  actually  corresponds  to  a  quadratic  function  rather  than  a  quartic.  This  fact  coupled 
with  the  use  of  2™^  order  finite-difference  approximations  provides  a  situation  in  which  exact 
derivatives  are  produced  on  an  equally  spaced  mesh.  Moreover,  it  is  observed  that  when  fc  of 
the  function,  f(x),  exceeds  that  of  the  parabolic  function,  points  tend  to  cluster  toward  the 
center  while  if  f^  is  less  than  that  of  the  parabolic  function,  points  cluster  toward  die  bound¬ 
aries.  This  information  may  be  used  in  future  mesh  adaption  schemes  to  improve  the  ability 
of  the  scheme  to  generate  a  “better  mesh"  with  assurance  of  improving  accurate. 

3.  Other  Pbssibilities 

In  theory,  one  could  determine  the  “best  mesh"  to  resolve  numerical  approximations  to  other 
types  of  terms  also,  such  as  the  quantity:  fx-uu„  found  in  the  viscous  Burgers’  equation.  In  this 
case,  however,  other  conditions  besides  no  inflection  points  must  likely  be  imposed  to  guar¬ 
antee  uniqueness  of  the  result,  if  such  uniqueness  is  even  possible. 

The  practical  usefulness  of  the  procedure  outlined  here  is  still  unknown  for  multi-dimension- 
al  multi-variable  problems.  However,  the  lessons  learned  are  clearly  valuable  guides  in  help¬ 
ing  to  formulate  new  weight  functions  and  new  adaption  schemes  which  incorporate  informa- 


38 


tion  regarding  both  the  nature  of  the  function  f(x)  and  the  discrete  approximation  used  to 
represent  its  derivatives. 


C .  2-D  Considerations 

New  considerations  are  important  in  extending  to  multi-dimensional  problems.  One  of  the 
most  obvious  is  the  issue  of  orthogonality.  The  basic  concern  is  the  error  introduced  into  a 
computation  as  a  result  of  non-orthogonality  of  the  mesh.  In  a  dynamically  adapting  mesh, 
this  concern  is  clearly  greater  since  the  real  possibility  exists  that  cell  collapse  may  be  inadver¬ 
tently  driven  by  the  adaption  scheme.  An  attempt  was  made  to  study  the  way  in  w^ich  non¬ 
orthogonality  influences  the  numerical  approximation.  The  basic  procedure  is  outlined  with 
the  key  results  highlighted. 

/  .  Errors  Associated  With  Mesh  Non-Orthogonality 


Consider  an  equation  of  the  form  u,  +  f,  +  gy  =  0.  An  application  ofa  generalized  mapping 
from  (t,x,y)  to  results  in  a  typical  term,  say  f*,  being  replaced  by  its  analytic  equivalent 

given  by;  f,  =  +  7*^7  or  f,  =  where  J  =  x^,,-yrpt^  .  If  the  quanti¬ 

ties,  f^  and  f,, ,  are  approximated  by  central  differences,  the  following  expression  results: 

ITie  next  step  is  to  apply  a  2-D  Taylor  expansion  to  each  f  term  and  then  collect  all  like  deriva¬ 
tive  terms  together  resulting  in: 


^7  y,;  (x%  1  -  Xj-i)  -  y|(xj  + 1  -  Xj_i)  |  ^ 

2JA|A>;  P 

A7  Vr,  [y,+  !-yi-i)-A^  yi(y,+  i-yy.i)] 

- - ^Vf« 

2JA^A7  ^ 

k. 


A7  y-7  (^4- 1  -  ^-i)  -  y|(^+ 1  -  ^.i)l 

4JA^A7  ’ 


^7  y?  I  -  yf-i)  -  ^  i  -  i^-i) 
4JA|A7 
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Let  this  be  ejqjressed  as  Cif*  +  Ciiy  +  C3fxx  +  C4fyy  +  Csf^y  +  ...  Then  it  is  desirable  to  have 
Cl  =  1  while  the  other  C’s  are  as  small  as  possible.  It  is  easy  to  show  that  Ci  =  landC2  =  Oif 
are  determined  from  the  central  difference  approximations: 

~  Xj+i-xj-i  ^  ^ 

^  ~  2A^  2A|  ’  ''  2^r|  2^r}  ’ 

~  ^'(y)  =  „  _  <^j(y)  _  yi^i-yn 

2A|  2A^  '  ~  2A7  2At) 

With  these  metric  evaluations  and  the  second  difference  definitions, 

(5f(x)  =  Xi+ 1  +  Xi-i ,  (52(x)  =  Xj  + 1  +  Xj_i ,  (5f(y)  =  y.+ 1  +  yi-i .  and  (5j^(y)  =  yj+ 1  +  yj-i , 
the  C3  and  C4  are  expressed  as: 

c  -  ~  <^»(y)^j(x)<^j^(x)  ^  ^  -  <?i(y)<?)(y)<3f(y) 

2[(5,(x)<5j(y)  -  <5/x)(5,(y)]  2[(5,(x)<5,(y)  -  (5j(x)<5,(y)] 

The  values  of  these  C’s  are  small  if  the  numerators  are  small  and/or  the  denominators  are 
la'ge.  Consider  first  the  C3  coefficient.  The  numerator  of  this  coefficient  is  zero  when  the 
equality,  p^  =  p,^  is  realized.  Inthisequality,  each  of  the  p’s  represents  a  percent  difference 

m  the  slope  of  the  corresponding  coordinate  line  between  the  +  and -cells.  For  example,  p^ 
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represents  the  difference  between  the  slope  of  the  rj  =  constant  line  connecting  i  with  i  +  1  and 
i  with  i-1  divided  by  the  slope  of  the  line  connecting  i  +  1  with  i-1.  This  is  expressed  as 

p  _  — hui  See  the  sketch  below. 

Si 


Since  this  is  not  likely  to  happen  for  any  mesh,  let  alone  a  dynamically  adaptive  one,  it  is  safe 
to  assume  that  Cj  will  be  non-zero.  Thus  it  is  even  more  important  that  the  denominator  is 
not  too  small.  Close  examination  of  the  denominator  reveals  that  it  is  proportional  to  the 
cross  product  of  the  two  vectors  connecting  i+  1  with  i-1  and  j+  I  with  j-1.  This  means  that 
the  denominator  is  proportional  to  the  product  of  the  magnitudes  of  these  two  vectors  and  the 
sine  of  the  angle  between  them.  This  is  where  non-orthogonality  can  be  detrimental  to  the 
accuracy  of  the  approximation.  As  the  mesh  tends  to  extreme  non-orthogonality,  the  denom¬ 
inator  tends  to  zero. 

Continuing  with  the  C4  coefficient,  it  is  clear  that  if  the  numerator  is  not  zero,  the  significance 
of  the  denominator  is  identical  tn  that  indicated  for  C3.  The  numerator  is  seen  to  be  zero  for 

the  most  general  case  when  df(y)  =  ^r(y).  Again,  this  is  not  at  all  likely  and  consequently  the 

effect  of  non-orthogonality  can  be  significant. 

The  Cs  coeffident  has  a  little  different  interpretation.  Since  it  multiplies  a  cross  derivative, 
one  would  expect  the  coefficient  to  have  .something  to  do  with  areas.  Cs  can  be  written  as: 

<^j(y)[j<,+  iyi+  1  -  x,_,y,,i)-  (5,(y)[:^+  iTj  +  l  -  ^iTj-i] 

[(5,(x)(5j(y)  -  <5/x)(5,(y)] 

From  the  following  two  sketches,  the  areas  Ai  -  A4  are  defined.  C5  is  then  expressed  in  terms 
of  these  areas. 
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C5  =  ~  ^3]  -  (3i(y){  A2  -  A4] 

[d,(x)<5j(y)  -  <5^(x)(5,(y)j 

From  this  interpretation,  it  is  clear  that  the  numerator  of  the  C5  coeffident  is  not  likely  to 
vanish,  except  in  special  situations,  and  consequently  the  significance  of  the  size  of  the  de¬ 
nominator  is  amplified.  Again,  this  depends  on  the  extent  of  the  non-orthogonality  of  the 
mesh. 

In  summary,  the  connection  between  non-orthogonality  and  error  has  been  established  for 
the  case  of  central  differences.  The  basic  results  are  expected  to  be  applicable  to  upwind  dif¬ 
ference  approximations  as  well  since  they  appear  to  be  common  sense  in  nature.  The  basic 
requirements  are: 

1.  metrics  should  be  evaluated  the  same  as 

2.  minimizing  coordinate  line  curvature  and  stretching  also  minimizes  numerators  of  Cj,  C4,  C5. 

3.  minimizing  non-orthogonality  maximizes  denominators  of  C3,  C4,  C5. 
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D  .  Dynamic  Adaption  For  2-D  Euler  Equations 

The  basic  grid  and  grid-speed  laws  for  2-D  problems  introduced  by  [  1  j  are  dimensionally 
inconsistent  unless  dimensions  are  absorbed  into  the  user  input  parameters.  The  cure  for  this 
scaling  problem  was  to  redefine  the  performance  index  in  a  dimensionally  consistent  manner 
with  user  input  parameters  of  ord*”’  1.  The  new  performance  index  is  given  as: 

I  =  1 1  L(r.  f^.  r,)d|d>;  ( 16 ) 

where  L  =  Ls  +  +  Lg  and: 


,  _  1  rt  •  Ff  +  f ,  •  r^ 

*-5  ““  j 


i-O  ~ 


2 


La  —  ^ 


K  = 

(^max  “  ^minX^nioX  “  7min) 

The  Euler-Lagrange  equations  were  then  re-derived  for  this  case  and  the  grid  speed  equa¬ 
tions  were  derived  by  differentiating  these  new  Euler-Lagrange  equations  with  respect  to  the 
temporal  variable. 

1 .  Modified  Mesh  Law 

The  result  of  the  derivation  was  a  second  order  PDE  system  of  the  form: 

G(r)  =  -l-  Cr^  +  Dr^  +  Er,  =  0  ( 17 ) 

which  constitutes  the  grid  law.  The  boundary  conditions  are  mixed  Dirichlet  and  Neumann 
depending  upon  the  phy  cal  conditions  required  of  the  mesh  at  the  boundaries. 

In  this  equation,  A=  As  -f  Aq  +  Ag,  B  =  Bj  +  Bo  +  Ba,  etc.  These  matrices  are  given  as 
follows: 
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This  new  mesh  law  requires  specification  of  3  parameters,  X*,  and  Xg  corresponding  to 
smoothness,  orthogonality,  and  adaption  terms.  Note  that  in  the  new  formulation,  these  may 
be  functions  of  time.  This  may  be  useful  for  problems  in  which  the  flow  structure  being 
adapted  to  diminishes  with  time  while  the  user  wishes  to  maintain  resolution  at  this  location 
during  the  decay  process.  These  parameters  are  nondimensional  and  have  the  value  of  unity  if 
orthogonality  and  solution  adaption  are  equal  in  importance  to  smoothness.  In  the  old 
scheme,  the  values  of  these  parameters  were  functions  of  the  problem,  the  number  of  mesh 
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points,  the  distribution  of  points,  etc.,  which  required  experimentation  with  various  values  in 
order  to  discover  an  acceptable  set  (if  one  even  existed). 

2 .  Corresponding  Mesh  Speed  Law 

The  result  of  the  temporal  differentiation  was  a  linear,  variable  coefficient,  second  order 
PDE  system  of  the  form: 


GT(f)  =  +  Cirjr,  +  D'zit  +  E'z,,  +  f  *  =  0  ( 18 } 

which  constitutes  the  grid  speed  law.  In  this  equation,  z  =  r, ,  the  grid  speed  vector.  The 
boundary  conditions  are  mixed  Dirichlet  and  Neumann  for  this  equation  as  well,  depending 
upon  the  physical  conditions  required  of  the  mesh  speed  at  the  boundaries. 

The  matrices  D*  and  E*  and  the  vector  T*  are  quite  complicated  and  involve  several  pages  of 
definition.  They  result  from  differentiating  A,  B,  C,  D,  E,  J,  Of.  3.  "V.  and  w  with  respect  to  t, 
and  then  expressing  the  result  as  a  matrix  times  the  grid  ^ed  vector  derivatives  where  possi¬ 
ble. 
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E .  Appendix  Figures 


Figure  (1).  f(x)  and  mesh  position  for  uniform,  “best”  central  and  “best”  one-sided 


Figure  (4).  f(x),  fx(x),  f„(x).  and  curvature  for  quartic(p  = -8). 
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Figure  (5).  f(x)  and  mesh  position  for  uniform,  “best”  central  and  “best”  one-sii 


Figure  (6).  f(x),  fj,(x),  fxx(x),  and  curvature  for  quartic(p  =  -4). 


